Crystalline phases of polydisperse spheres 
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We use specialized Monte Carlo simulation methods and moment free energy calculations to 
provide conclusive evidence that dense polydisperse spheres at equilibrium demix into coexisting 
fee phases, with more phases appearing as the spread of diameters increases. We manage to track 
up to four coexisting phases. Each of these is fractionated: it contains a narrower distribution of 
particle sizes than is present in the system overall. We also demonstrate that, surprisingly, demixing 
transitions can be nearly continuous, accompanied by fluctuations in local particle size correlated 
over many lattice spacings. 



Suspensions of spherical colloids have long served as 
an experimentally accessible testing ground for our un- 
derstanding of the liquid, crystalline and glassy states 
of matter [1, 2]. Such work is complemented by theory 
and simulation, which attempt to reproduce, rationalize 
and predict experimental results. In so doing, it is com- 
mon to treat the suspension as an assembly of identical 
spheres. But this neglects a key feature, namely that 
the chemical processes by which real colloids are synthe- 
sized invariably produce particles that have a spread of 
diameters, i.e. they are 'polydisperse'. As is becoming 
increasingly clear, polydispcrsity gives rise to a rich va- 
riety of novel phenomena not observed in monodisperse 
systems [3]. However, despite sustained attention, basic 
questions remain concerning its effects on one of the most 
fundamental aspects of any thermal system, namely the 
equilibrium phase behaviour. 

A case in point is the character of the thermodynam- 
ically stable structures of size-disperse spheres in the 
dense regime, above typical fluid densities. Polydisper- 
sity should act to destabilize a crystal because of the dif- 
ficulty of accommodating a range of particle sizes within 
a single lattice structure; but there has been no definite 
answer as to what stable structures arise instead. Indeed, 
the nature and extent of the influence of polydispersity 
both on the crystalline phases and the location of the 
freezing line is controversial. On the theoretical front, 
there is a diverse range of predictions of novel phenom- 
ena including reentrant melting [4], an 'equilibrium glass' 
phase [5], and solid-solid coexistence [6-8]. Addition- 
ally, recent simulation work has reported the occurrence 
of a partly crystalline 'inhomogeneous phase' within an 
approximate phase diagram based only on equality of 
single-phase free energies [9, 10]. Other simulations sug- 
gest that the fluid-solid coexistence region terminates in 
a critical point beyond which a disordered solid occurs 
[11]. On the experimental side, studies of colloidal sys- 
tems observe that beyond a certain 'terminal' polydisper- 
sity no crystallization occurs on experimental timescales 
[1], although it remains unclear whether this is a true 



equilibrium effect or a manifestation of dynamic arrest. 

A crucial distinction between monodisperse and poly- 
disperse systems at phase coexistence is the ability of the 
latter to fractionate so that the distribution of the par- 
ticle diameters, a, varies from phase to phase [12-14]. 
If for a certain phase (labeled a), one counts the num- 
ber density of particles having diameters in the range 
a . . . a + da, this serves to define a density distribution 
p( a \a). Experimentally, however, for most complex flu- 
ids one has the constraint that the overall distribution of 
sizes (across all phases) has a form fixed by the synthe- 
sis of the fluid. This gives rise to a generalized lever 
rule: p^°\cr) = J2a^ a P^( a )> w ^ n ^ a the fractional 
volume occupied by phase a, p(°\a) the 'parent' den- 
sity distribution and p^ (a) the 'daughter' distributions. 
Since the form of the parent is fixed, only its scale is 
free to vary, e.g. by dilution with solvent, and one writes 
p(°\a) = n°f(a), where n° is the total number density 
and f(a) is a prescribed normalized shape function. The 
polydispersity, S, is then defined as the standard devia- 
tion of the parent distribution, in units of its mean. 

The diversity of theoretical and simulation findings 
stems from the sensitivity of the results to the accuracy 
with which fractionation is treated. Previous work has 
either disregarded fractionation entirely, or used drastic 
(and differing) approximations to describe it. An excep- 
tion are moment free energy (MFE) theory calculations, 
which do account fully for fractionation, and which have 
previously been reported by one of us for hard spheres 
[15]. These predict that increasing polydispersity shifts 
the fluid-solid coexistence region to higher number den- 
sities, but that neither reentrant melting nor a terminal 
polydispersity occurs. Instead, the fluid can always split 
off a small volume of dense phase whose size distribu- 
tion is sufficiently narrow for crystallization. Moreover, 
as one increases n° or S within the solid region, a suc- 
cession of phase transitions is predicted in which the sys- 
tem demixes into an ever greater number of differently- 
fractionated 'daughter' phases. However, the MFE cal- 
culation uses approximate free energy expressions, which 
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for solids are derived from those of binary mixtures and 
implicitly already assume that all solids are fee. Inde- 
pendent confirmation of its predictions is then highly de- 
sirable, but has hitherto been lacking. In this Letter we 
provide a definite answer to the question of the nature 
of the equilibrium phase behaviour via state-of-the-art 
Monte Carlo (MC) simulations, and compare with MFE 
calculations; both fully provide for fractionation and em- 
ploy a fixed parent size distribution. 

In simulations the appropriate framework for observ- 
ing genuine equilibrium behaviour in dense polydisperse 
particles is the isobaric semi-grand canonical ensemble 
[16, 17]. This is the analog of a monodisperse (N,p,T) 
ensemble where the prevalence of different particle sizes 
is controlled by imposing in addition chemical poten- 
tial differences A/x(<r) that are measured relative to the 
chemical potential of some reference particle size. Monte 
Carlo sampling of this ensemble can exploit particle re- 
sizing moves to allow local sampling of the size distribu- 
tion without the need for particle diffusion (thus catering 
for fractionation effects) , while volume updates facilitate 
density fluctuations so that the system can transform be- 
tween phases. Our study is the first to deploy this ensem- 
ble in the crystalline regime together with a method for 
imposing a fixed overall parent distribution. This allows 
determination of physically realistic phase behaviour in- 
cluding the boundaries of the onset of coexistence (known 
as cloud curves) and daughter distributions. Additionally 
we can calculate - but do not show here - shadow curves 
which record the density and volume fraction of the new 
phase when coexistence first occurs. Cloud and shadow 
curves do not coincide, demonstrating further the pres- 
ence of fractionation: new phases that appear generically 
have size distributions different from the parent [3]. We 
combine the above techniques with the specialized phase 
switch Monte Carlo (PSMC) method [18, 19] for obtain- 
ing fluid-solid coexistence properties. In both cases, the 
chemical potential differences Afi(a) are determined it- 
eratively to match the ensemble-averaged density distri- 
bution (p(u)) to the prescribed parent p^\a) = n°f(a). 
At coexistence, this is supplemented by an equal peak 
weight criterion for the order parameter distribution to 
ensure that finite-size effects are exponentially small in 
system size [20, 21]. 

We stress that the choice of ensemble and use of sophis- 
ticated sampling and analysis techniques are crucial to 
observing qualitatively correct phase behaviour in poly- 
disperse systems. Use of standard canonical [9] or mi- 
crocanonical ensembles [10, 11] are unequal to the task 
and almost certainly yield major artifacts. The reasons 
for this are three fold: (i) the dynamics is too slow to al- 
low fractionation on simulation timescales; (ii) the sizes 
of the particles are fixed, which for a finite system pre- 
vents daughter distributions assuming an arbitrary form 
as they can in the thermodynamic limit; (iii) these ensem- 
bles necessarily form interfaces between coexisting phases 



and for accessible particle numbers one cannot hope to 
see multiple coexisting crystalline phases when this oc- 
curs. 

Our simulations consider a system of 256 particles in- 
teracting via a strongly repulsive pair potential 

v{r ij )=e((Ti j /ri j ) 12 , (1) 

with particle distances = \n— rj \ and interaction radii 
Cij = (<Ti+aj) /2. The choice of this potential rather than 
infinitely repulsive (hard) spheres is made on pragmatic 
grounds: an MC contraction of the simulation box that 
leads to an infinitesimal overlap of two hard spheres will 
always be rejected, so (particularly at high densities) we 
can expect higher MC acceptance rates using this 'softer' 
potential. In common with hard spheres, the monodis- 
perse version of our model freezes into an fee crystalline 
structure [19, 22], and temperature only plays the role of 
a scale: the thermodynamic state depends not on n° and 
T separately but only on the combination n (e//csT) 1 / 4 . 
Phase diagrams for different T then scale exactly onto 
one another, and we can fix e/ksT = 1. 

We consider parent size distributions of the top-hat 
form: 

f(a) = { (2c)- 1 ifl-c<a<l + c 

■' x ' [0 otherwise v ; 

The width parameter c controls the polydispersity 5 — 
c/-\/3, and we have set the mean particle diameter to 1. 
With these choices, and the interaction potential (1), our 
results are directly comparable to the phase diagram of 
Ref. [9] where neither fractionation nor, at a more basic 
level, the presence of coexistence regions of finite width 
was allowed for. 

Using PSMC we have mapped the cloud curves of the 
fluid-solid transition, using S as our control parameter, 
up to polydispersities of 6 w 8.7% on the fluid side and 
8 ss 7% for the solid, both of which are in the typical 
range for colloidal systems, but not so great that we ex- 
pect to see exotic phases such as ABi 3 . As shown in 
Fig. la, both the fluid (circles) and the solid (squares) 
phase cloud densities shift to higher n° as S is increased, 
but without the sharp narrowing that would be required 
for a reentrant melting scenario [4] . 

Turning now to the solid region, a comprehensive ex- 
ploration of the (n°-S) plane is impractical because of 
the relatively high computational cost of our specialized 
simulation technique. But we can understand important 
qualitative features by following the dashed trajectory 
included in Fig. la. Along this path, we monitored the 
state of the system via the probability distribution of the 
fluctuating total number density p(n), which serves as 
an order parameter for phase changes. Starting from the 
fee solid cloud point at 6 = 6.3%, we initially increased 
n° in a stepwise fashion (filled circles) to n° = 1.45, and 
then switched to increasing 5 at constant n° as a po- 
tentially faster route to demixing. Indeed, at 5 w 8% 
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FIG. 1: (Color online), (a) Simulation results for the partial 
phase diagram of the model (1) with parent distribution (2). 
Asterisks: points where new solid phases appear; dashed 
lines: phase boundary slopes found by histogram reweighting. 
F=fluid, S=solid. Colored symbols: state points considered 
in Fig. 2. (b) MFE calculation of phase diagram of hard 
spheres with the same parent form. The dashed line shows a 
trajectory comparable to that followed by the simulations. 



there was a smooth change in p(n) from single to double 
peaked; an example of the double peaked form is shown 
in Fig. 2. The two associated phases were identified as 
being fee solids. As is physically reasonable, the higher 
density solid (HDS) daughter phase contains a surplus of 
the smaller particles while the lower density solid (LDS) 
phase has more of the larger particles. 

Continuing to higher S eventually led to spontaneous 
melting of the system at S = 13.7%, implying that the 
limit of metastability with respect to a fluid-solid-solid 
(FSS) coexistence had been overstepped, as is indeed pre- 
dicted by our MFE calculations (see below). We there- 
fore backtracked slightly into the solid-solid (SS) region, 
embarking on a new trajectory with increasing n° at con- 
stant S — 13.5%. This produced a third peak in p{n) 
at n° ss 1.475. The corresponding intermediate density 
solid (IDS) was again found to be isostructural with the 
other two, with dominant particles sizes between those in 
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FIG. 2: (Color Online). Distribution of the overall number 
density p{n) at the SS and SSSS statepoints indicated by the 
colored symbols in Fig. la. 




FIG. 3: (a) Solid: parent distribution at (n° = 1.73,5 = 
13.5%). Symbols: Simulation results for the four daughter 
distributions. The associated fractional volumes \ a are (left 
to right) 0.209, 0.188, 0.232, 0.373. (b) MFE results at the 
comparable state point (n° = 1.232, 8 — 8.7%); fractional 
volumes are 0.273, 0.162, 0.200, 0.365. 



the HDS and LDS. Finally, increasing the overall density 
to n° ~ 1.68 we observed that the central IDS peak in 
pin) split rather smoothly into two peaks, yielding a four 
peaked structure (Fig. 2): four fee solids now divide the 
range of particles sizes among themselves (Fig. 3). 

We next compare to our theoretical MFE calculations. 
These used the same parent size distribution (2) but, 
since no suitable polydisperse model free energies are 
available for the soft repulsive potential (1), the analysis 
was performed for hard spheres, using the methodology 
described elsewhere [15]. The qualitative physics should 
be the same. Indeed, taking a comparable path [23] 
through the calculated phase diagram (Fig. lb) shows the 
same features as in the simulations. (Quantitatively, the 
fluid-solid coexistence region is narrower, and transitions 
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to multiple solids occur at lower n° and S, presumably be- 
cause with a hard repulsion, a crystal can accommodate 
above average-sized particles less easily.) Also the frac- 
tionation effects are well reproduced, as shown in Fig. 3b 
for an SSSS state point at a location comparable (relative 
to phase boundaries [23]) to the one in Fig. 3a. 
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Our tailored simulations have provided a clear an- 
swer to long standing questions surrounding the effect 
of size polydispersity on the equilibrium phase behaviour 
of spherical particles: as density and/or polydispersity 
are increased within the crystalline region, the system 
demixes into an ever increasing number of fractionated 
fee phases. Given the high level of qualitative accord 
with MFE calculations, we are confident that this sce- 
nario represents the true equilibrium situation. Since 
the MFE results are insensitive to whether the parent 
distribution is top hat (present work) or has a Schultz or 
triangular form [15], we believe the demixing scenario to 
be quite general. Understanding in detail when and why 
the clemixing transitions are near-continuous is an excit- 
ing open challenge. Finally, we note that in spherical 
colloids demixing transitions may not always be directly 
observable because fractionation requires particle diffu- 
sion which is inhibited in solids. Nonetheless, one might 
expect to see evidence for solid demixing in regions where 
the solids coexist with a fluid (cf. Fig. lb) that can trans- 
port particles to their preferred solid phase. Even in sit- 
uations where equilibrium cannot be reached for kinetic 
reasons, knowledge of the true equilibrium state provides 
an important baseline for interpreting dynamical effects 
[2, 24]. 



FIG. 4: Correlation volume £ 3 in the solid phases encountered 
along the phase diagram trajectories of Fig. 1. (a) Simula- 
tions, (b) MFE calculations. 

A surprising feature of our results is that, from the vari- 
ation of p(n), the transitions S— >SS and SSS^SSSS ap- 
pear to be near-continuous in character, while SS— t-SSS is 
strongly first order as is usually expected for transitions 
in the solid state. A near- continuous transition should 
be accompanied by size fluctuations correlated over large 
distances, as precursors of the new phases, whereas the 
fluctuations will remain small on approaching a first or- 
der transition. To quantify these fluctuations we measure 
a correlation volume £ 3 from the variance across configu- 
rations of the mean particle size a. Suitably normalized, 
this variance, ((Act) 2 ), is proportional to the spatial in- 
tegral over the pair correlation g aiJ >{r), weighted by de- 
viations of the particle sizes a and a' from the mean, 
i.e. the correlation volume. In theoretical calculations, 
((Act) 2 ) can be extracted from second derivatives of the 
MFE [13]. Measurements of £ 3 along the trajectories 
through the phase diagrams are shown in Fig. 4. This 
grows large near the transitions to two and four solids, 
confirming their near- continuous character. In the latter 
case, the splitting of the middle peak seen earlier in p(n) 
suggests that the new solids arise out of the IDS phase, 
and this is consistent with large fluctuations occurring 
(see Fig. 4) only in this phase and not the HDS or LDS. 
The MFE predictions are, again, in good qualitative ac- 
cord with the simulation data. 
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